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Collapse, or a gravitational-like phase transition is studied in a microcanonical ensemble of 
particles with an attractive 1/r" potential. A mean field continuous integral equation is used to 
determine a saddle-point density profile that extremizes the entropy functional. For all < a < 3, a 
critical energy is determined below which the entropy of the system exhibits a discontinuous jump. 
If an effective short-range cutoff is applied, the entropy jump is finite; if not, the entropy diverges 
to -l-cxD. A stable integral equation solution represents a state with maximal entropy; the reverse is 
always true only for a modified integral equation introduced here. PACS numbers: 02.30.Rz 04.40.-b 

05.70.Fh 64.60.-i 



I. INTRODUCTION 

The behavior of systems with long-range interactions is often different from those considered in traditional ther- 
modynamics. As an example we take self gravitating systems, i.e. particle ensembles with pm'ely attractive 1/r 
interactions. It is known that such systems exhibit collapse, sometimes called a zero-order phase transition, when the 
energy in the microcanonical ensemble (ME) or the temperature in the canonical ensemble (CE) drop below a certain 
critical value Cc or Tc, respectively During such a transition, the corresponding thermodynamic potentials 

(entropy in the ME or free energy in the CE) exhibit a discontinuous jump. If the interaction between the particles is 
purely attractive and no short-range cutoff is introduced, then the discontinuous jump is infinite and the entropy and 
free energy go to +00 and —00, respectively. This makes all normal (noncollapsed) states of the self-attractive system 
metastable with respect to such a collapse; the collapse energy Cc is in fact an energy below which the metastable 
state ceases to exist. If, on the other hand, some form of short-range cutoff is introduced, the entropy and free energy 
jumps are finite. In this case as a result of the collapse, the system goes into a nonsingular state with a dense core, 
the precise nature of which depends on the details of the short-range behavior of the potential. Then only the normal 
states which are in some interval of energies above the collapse point are metastable with respect to such a transition 
(see Fig. 1). 




e 

Fig. 1. Sketch of an entropy vs. energy plot for a system with gravitational- like collapse. The entropy of the normal (non- 
collapsed) state is shown by a solid line, the entropies of the two collapsed states for different cutoff radio ri and r2, r2 < ri, 
are shown by dashed lines. The entropies of the two collapsed states intersect the entropy of the normal state at energies 
and 62. 

There is an energy e* for which both collapsed and normal system have the same entropy (cf. e* and in Fig. 1); 
above this energy the collapsed state becomes metastable and at some higher energy it ceases to exist [^jj^. It is 
possible to regard the energy e* as that where a true phase transition occurs. When the effective cutoff vanishes, such 
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phase transition energy e* increases to infinity, so that without a cutoff all the finite energy states are metastable Q . 
However, the value of e* is highly sensitive to the details of the short-range cutoff. On the contrary, the collapse energy 
Cc depends on the long-range part of the interparticle interactions and is almost unaffected by a cutoff, provided that 
it is sufficiently short-range. Since the collapse phenomena carry most of the information about the peculiarities of the 
phase behavior of the system caused by the long-range interactions, they usually gain the most attention, considering 
then a phase transition at e* as a mere artifact of an often auxiliary cutoff. 

While these rather elaborate studies of gravitational collapse have usually been motivated by cosmological applica- 
tions and were performed solely for an 1/r potential, a natural question is, what happens if the particles interact via an 
attractive l/r" potential with arbitrary a. It has been noticed before that in systems with nonintegrable interactions, 
i.e. when a is less than the dimensionality of the space, first order phase transitions occur differently in the ME than 
in the CE, even for iV — > cxd However, in the examples considered in the literature, the potential energy was 
always bounded from below (usually by putting the system on a lattice), allowing only for normal first-order phase 
transitions and excluding any singular collapse. 

In this work we consider a possibility of collapse in systems, similar to the gravitational self- attracting Hamiltonian 
particle systems, but with a general 1/r" potential. The paper is organized as follows: after this brief introduction we 
formally define the model and derive an integral equation to solve for a saddle-point density profile. Then we describe 
the results obtained by a numerical iterative solution of this equation for various < a < 3. After this, we look 
at the local stability of the solutions of the integral equation and link their stability to the stability (metastability) 
of the corresponding states of the self-attracting system. Finally, we discuss the results and directions for further 
investigation. 

While the collapse is shown to be a generic property of self-gravitating systems occurring in both the ME and the 
CE, the ME allows to obtain more information about the system than the CE. For example, negative heat-capacity 
states, that precede the collapse in the ME, are not accessible in the CE, as it is in the much better studied case 
of normal first-order phase transitions |^,0j9j. In the rest of our paper we consider systems only in the ME. In our 
studies we strongly rely on the large body of existing results derived for a = 1, often rigorously p|,p|JTc[|. 



II. THE MODEL 



Let us consider a 3-dimensional ME of N particles each having mass m, interacting with an attractive pair potential 
Vij = —G/\ri~rj\°-. These particles are confined to a spherical container with radius i?; the total energy of the particles 
is E. The microcanonical entropy (equal to the logarithm of the density of states; here and in the following we put 
ks = 1) can be expressed as: 
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After integrating over the momenta and introducing a dimensionless energy e = ER!^ /GN"^ and radial coordinates 
Xi = Ti/R, Eq. (lit) can be reduced to: 



S{e) = 5-0 + log 
where Ckm is a dimensionless kinetic energy. 
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and 5*0 is an energy-independent term. 
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Here T(x) is a Gamma function, Q(x) is a unit step function which guarantees positiveness of the kinetic energy, 
and the integration for each dxi runs over a 3-dimensional unit sphere. It follows from Eq. (|^) that the entire 
thermodynamical behavior of our system depends on the single variable e; and for a system to remain in the same 
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state when the number of particles is varied but the energy per particle is fixed, the system size R should scale as 
R ~ N^/". This indicates that a = 3 is the largest value of a for which no thermodymanic scaling is observed and 
hence some non-traditional phase transitions may exist. 

We assume that the number of particles N is large and go to a continuum limit, replacing a 3Af-dimensional 
configurational integral in (0) by a functional integral over possible density profiles p{x) (see, for example, 0,01): 



5(e) «5o + log |y Vp 
where the "effective action" functional s is defined 
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Here we introduced two auxiliary Fourier integrations: one over d/3 to replace the Q function in Eq. (0) (see e.g. 
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and similarly one over d'y to express the S function S [J p{x)dx — l] . The Eq. (^ is applicable only for < a < 3 
when the integrals on dxi and dx2 are convergent at the lower limit. 

Using that N is large, the integral in Eq. (|^) can be evaluated by the saddle-point method. To determine the 
dominant contributions to the integral, we differentiate (^ with respect to 7, /?, and p{x) and look for {7s, /3s, /5s(x)} 
which give the extrema to the effective action (^). As a result, unique solutions are obtained for {7s,/?s}; while for 
Ps{x) a nonlinear integral equation is obtained for which the number and the nature of solutions is generally unknown: 
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It follows from (||) that is equal to 3/2 of the kinetic energy, hence it is equal to the entropy derivative with respect 
to the energy, i.e. to the inverse temperature, Ps — ds{e)/de = l/T. 

Using (^Q and ignoring terms independent of e, we obtain for the entropy per particle s(e) = S{e)/N in the 
saddle-point approximation: 



In 



p{xi)p{x2] 

\X1 - f2|" 



dxidx2 



p{x)\n[p{x)]dx + 0{^). 



(9) 



The saddle-point approximation works only if the second derivatives of the effective action functional s are not too 
small; we will return to the question of its validity in Section IV. 

The above equations look similar to those for the gravitational (a — 1) case which are derived, for example, in [^,^, 
to where the reader is referred for a more detailed description of the derivation. However, we can not proceed further 
along the lines developed in the references mentioned above, which all deal with the gravitational case. For, unlike 
the gravitational case, when the A(l/r) — — 47r(5(r) property of a Coulomb interaction 1/r allows one to reduce (^) 
to a second-order differential equation, no such property exists for a general a, and we have to deal with the integral 
equation for Ps{x) in (||) directly. 

In the following a further simplification is made by disregarding the angular dependence for the saddle-point density 
profile, i.e. p{x) = p{x)] for the gravitational case this is justified in 1^,^. After performing a straightforward angular 
integration in the exponential of Eq. (@), we obtain for a ^ 2, 
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and for a = 2, 



exp 



(2 - a)x' Jo 
1 /.I 



dx' 



dxidx2 



= po exp ^ / p{xi)xi^n{x + Xi) — \n\x ~ xi\\dxi 

X Jo 



Po 





■27r/3, /■ 


[Jo 


[ ^' Jo 







p{xi)xi (ln(a;' + xi) — In \x' — Xi\) dxi 



dx' 



Jo 



p{xi)p{x2)xiX2 [ln(.T + xi) — In \x — xi \] dxidx^ 



(10) 



(11) 



Once the Eqs. (10, Tl|) have been solved and the entropy (^ has been calculated, the nature of the phase behavior of 
the system can be deduced from an entropy-energy plot. In the next section we analyze the Eqs. ( [l0| ^ [Tl] ) numerically. 



III. NUMERICAL SOLUTION 

If /3s is considered an independent parameter rather than a factor depending on e and ps, the Eq. (||) is often 
called a Generalized Poisson-Boltzmann-Emden Equation ||l^ . Very little is known about this equation even in the 
gravitational a = 1 case. The only exactly known, so called "singular", solution is for /^^ = 2 and has the form 
Psing{x) = {Attx^)'^ , which leads to e = —1/4 and s = ln(\/27r) — 2. However, we are interested here in general a and 
e; and Ps in the ME is not an independent variable but a function of Cs- As we will discuss below, this has a crucial 
effect on the numerical accessibility of solutions in a certain interval of e. 

To solve the Eqs. ( p^ , pT| ) numerically for general a and e, we use a simple iterative method. For a fixed a we start 
at a relatively high value of e with a flat density profile, po = 3/47r, which is the solution of (|l^, |ri|) for e +00. 
Putting Po into (|lO|, 111]), we calculate /3s and obtain a new density profile p\(x). In other words, an iterative map 



Pi+\(x) = F^{p^{),x\ 

is introduced, with a nonlinear functional Fe[pi(.)j_x] defined by (|l^, [Tl| ). 
After a sufficient convergence of the iterations ( |l2| ) is achieved, i.e., when 
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the entropy is calculated with (^. We then move to a lower energy point, use the previous energy point density 
profile as poi repeat the procedure again. Caution in selecting the initial energy as well as the energy step has 
to be exercised to maintain the positiveness of the kinetic energy. In fact the initial energy has to be larger than 
~ (6~aK4~a)(3-a) ' '^l^i^h is the potential energy for the flat density profile. As we progress towards lower energies, the 
density profiles get more and more peaked near the origin, and the absolute value of the potential energy increases. 
To allow for improper integrals in the numerical integration of (|l^, pi] ), we use a simple mid-point trapezoid rule. 
Uniform meshes of 1000-2000 points were usually sufficient. However, in order to determine the position of a phase 
transition point with a sufficiently high precision in order to compare our results to the existing ones for a = 1, and 
to achieve also a discontinuous phase transition for smaller a, we had to use finer meshes with up to 5000 points. The 
convergence parameter 5 was usually set to be 10~^. 

The main conclusion that can be derived from the numerical results obtained is the following: for all < a < 3, as 
for a = 1, there is a certain energy tc{ot) below which the system collapses and the entropy exhibits a discontinuous 
jump. The results for tc{ol) are presented in Fig. 2. 
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Fig. 2. Plot of collapse energy ec{a) vs. potential exponent a. 

To verify our calculations of ec(a), we compare our result for eda — 1) with the existing data obtained by other 
methods. Our number, ec{a = 1) = —0.3346, is consistent with eda = 1) = —0.335, quoted in |^j4|. 

Since the behavior of self-attractive systems is quantitatively similar for all < a < 3, let us consider in more detail, 
for example, a system with a = 1/2. Plots of entropy s(e) and inverse temperature /3s(e) = ds{e)/de are presented in 
Fig. 3. 
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Fig. 3. Plots of entropy s(e) and entropy derivative /3s(e)) for a — 1/2. The radius of excluded central volume ro = 5 x 10 
The points ec, e*, tu, and are defined in the text. 



As we go down along the energy axis e, the entropy decreases as well, passing through an inflection point 
where /3 reaches its maximum /J^. For energies below this inflection point, the system has a negative specific heat 
{cPs{e)/de'^ > 0) and is therefore unstable in the CE. As we pass through the point and continue decreasing the 
energy, the convergence of ( |T^ ) becomes slower and slower, and at the point Cc the iterations start to diverge. It is 
straightforward to show for all < a < 3 (see, e.g. for a = 1) that the entropy is unbounded from above with 
respect to uniform squeezing of all the matter into a sphere with a radius going to zero. Hence, if no short-range 
cutoff is present, it is reasonable to assume that the entropy discontinuity at eda) is infinite. 

If some form of a short-range cutoff is introduced, the entropy discontinuity may become finite. To investigate this 
we tried two approaches. One, suggested in is to place a small spherical excluded volume with a radius ro in the 
center of the system, or, in other words, to replace a spherical container with a spherical shell container. The other 
approach is to replace the original "bare" potential 1/r" with a "soft" potentials of the form l/(r^ -I- Tq)^". For a 
reasonably small short-range cutoff (rp ^ lO"'^ for small a, tq ~ 10^^ for a ~ 3 for both approaches) the behavior 
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of the pre-collapsed system is virtually unaffected. A finiteness of the integration mesh can also play the role of 
an effective central excluded volume; for our method of integration the size of such an effective excluded volume is 
roughly of the order of the mesh step. 

Introduction of a short-range cutoff makes the existence of a non-singular collapsed phase possible. However, being 
applied directly, the iterative method ( |T^ ) still diverges when e < ec- To make it convergent, we introduced a map 
with a variable "step" , 

p,+i{x) ^ <7F[p,i.),e] + {1 - a)p,ix), (14) 

where < tr < 1 is the step size parameter. Choosing a sufficiently small (as small as ^ 10^^ — 10^'^), we were able to 
make the algorithm convergent for e < Ec- The connection between the numerical stability of the iterative algorithm 
and the thermodynamic stability of the corresponding phase is analyzed in the next section. A typical density profile 
in the collapsed phase exhibits a much higher concentration around the origin than to the normal (un-coUapsed) 
phase; plots of density profiles for a = 1/2 are presented in Fig. 4. 
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Fig. 4. Density profiles p{x) for a = 1/2 for normal (dashed line) and collapsed (solid line) phases for the energy e* = —0.708, 
when entropies of both phases are the same. The radius of the excluded central volume is ro = 5 x 10~^. 

A collapsed phase exists not only for e < ec, but for e > ec as well. In fact, this phase is globally stable in the range 
of energies where its entropy is higher than that of the normal phase, i.e. when e < e*. For e > e*, the collapsed phase 
is metastable and above some energy e„ becomes unstable even locally (see Fig. 3). However, being cutoff-dependent, 
e* and e^ are not directly related to the fundamental properties of the original 1 /r" self-attracting system. 

Finally we return to the exact Psmgi^) = (47ra;^)^^ solution which exists for e = — 1/4 and a = 1. Our attempts to 
approach this solution by the numerical iterative methods ( ^^|Jl^ ) failed. In fact, even after substituting the psing{x) 
into ( [T^ ) as an initial approximation pq{x), the iterative solution of ( [T^ ) evolved either to a normal or to a collapsed 
solution depending on the value of the step a. We calculated the entropies for the three solutions, that exist at 
e = — 1/4: a normal s„, a collapsed Sc (with a sufficiently small central core, so that it still exists at this energy), and 
a Ssing for Psing- It turns out that Ssing < min{sn, Sc}, which, together with the evidence obtained from the iterative 
procedures mentioned above, strongly suggests that in the space of solutions (or fixed points) of (p^), both normal 
and collapsed p{x) are at least locally stable (attractive), while is unstable (repulsive). 




IV. STABILITY OF THE THERMODYNAMIC STATE AND THE ITERATIVE MAP 

The iterative solutions of the integral equation ( |l0| , |l2|) for the saddle point density profile ps {x) of the self-attracting 
system correspond to thermodynamically stable or unstable states. In this section we will give the necessary and the 
sufficient conditions for the thermodynamic stability of the Ps{x) in terms of the stability of the iterative solutions 
of the integral equation (p^). To that end we will look at the stability of ( p^ for a certain trial function pi{x). Let 
us assume that Ps{x) is a solution (fixed point) of (p^, and Spi{x) = pi{x)— Ps{x) is a small deviation of the ith 
iteration from the solution Ps{x). After one iteration we obtain for 5pi^i{x) = pi^i{x) — Psix) 
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For the iterations to converge, |5pi+i(a;)| < |5pi(x)| for all x. This condition is equivalent to a requirement that the 
absolute values of all the eigenvalues of the first variation operator defined below are less than one, 
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where the 4>i are the eigenfunctions corresponding to A^. Recalling the definition of F(:[p(.), x] (|8|, |lCl| , |ll| ) we note that 
it can be expressed through a first variation of the effective action functional s[p(.), e, 7, /?]: 
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Then one obtains for the first variation of Ff^[p{.)^x\: 
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Now for a state ps {x) to be thermodynamically stable and for the saddle point approximation to be applicable at all, 
the second variation of the effective action s on the right hand side of ( p^ must be negative for all x and x' . This is 
equivalent to the requirement that all the eigenvalues of the operator on the right-hand side of ( [l8|) are less than one. 
For, the only eigenvalue of the delta- function 5{x — x') is one with any function being its eigenfunction and the density 
Ps {x) is strictly positive for all < a; < 1 . Hence, the convergence of the map (|l^) to the function ps {x) is a sufficient 
condition for the thermodynamic stability of Ps{x) and the validity of the saddle point approximation. However, the 
reverse is not true, i.e. stable thermodynamic states do not necessarily correspond to stable iterative solutions. For, 
the operator on the left-hand side of Eq. ([l^ ) may have eigenvalues which are less than -1, which will make the map 
( |l^ ) unstable. Therefore let us instead of (j^ consider the variable step map (|lj) with eigenvalues of its first variation 
A*((t) = a\i -\- {1 — a). Then for any negative A^ we can chose an appropriate a so that |A*((t)| < 1. Evidently, this is 
what happened in the case of the collapsed phase which, possessing the highest entropy, is definitely stable, but can 
be accessed numerically only by using (|lj) with a sufficiently small a. On the other hand, for thermodynamically 
unstable states, such as Psing-, some of the eigenvalues of the operator in the right-hand side of Eq. (|l8|) are larger than 
1, which make the iterative maps (|ljjlj) unstable, so that such solutions cannot be found iteratively. This completes 
the demonstration of the equivalence of the two stabilities. 



V. DISCUSSION 



1. In the previous sections we revealed the existence of collapse and associated with it a discontinuity in the entropy 
in the microcanonical ensemble of particles with 1/r" attraction for < a < 3. This discontinuity was an infinite jump 
if no short-range cutoff was present. A carefully introduced short-range cutoff leaves the properties of pre-collapsed 
system virtually unaffected, but makes the entropy jump finite and allows to observe the collapsed phase. We proved 
that the stability of a solution of the integral equation for a saddle-point density profile is a sufficient condition for 
the profile to make the entropy a maximum and therefore be either stable or metastable. We modified the integral 
equation to make the reverse also true. 

2. In the range of a we have been working with, < a < 3, the potential is often called "non-integrable" , since the 
integral / d^r/r°' diverges at its upper limit. As the potential becomes integrable {a > 3), the continuum approach 
used here becomes inapplicable because the short-range density fluctuations, which the continuous approach cannot 
account for, become dominant over the long-range effects. Formally, the short-range nature of the behavior of the 
systems for a > 3 manifests itself as the divergence of the integral / d^r/r°' at its lower limit. 

3. A very important question which remains is that of the order of the gravitation- like phase transition. Here we 
have to distinguish between the collapse itself, which happens at Cc, and the phase transition which happens at the 
energy e* where the entropies of non-collapsed and collapsed states are equal (see Fig. 3). Since the entropy at the 
collapse point exhibits a discontinuous jump, the collapse is often called a zero-order phase transition However, 
the collapse is not a phase transition in the formal sense since it converts a metastable state into a stable one, which 
can be singular or finite, depending on the presence of the short-range cutoff. 
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On the other hand, the "true" phase transition between stable phases, which happens at e*, is sometimes referred 
to as a "gravitational first-order phase transition" j|] . Its distinct features include an inability of the two phases (non- 
collapsed and collapsed) to coexist and a discontinuous /3(e) i.e. temperature ||]. Yet in the "normal" ME first-order 
phase transition in a long-range interacting system (such as a mean- field Potts model), /3((e) remains continuous and 
smooth, but exhibits non-monotonous behavior: the interval of energies where phases coexist includes an interval 
where d(3{e)/d€ is positive and the specific heat is negative (see Fig. 5) [0-|9|. Hence there is an intrinsic difference 
between the normal and the gravitational first-order phase transitions. 

Even more, normal first-order phase transitions are found to replace gravitational first-order phase transitions which 
occur in the self-attracting systems that we consider here, if the short-range cutoff is sufficiently increased. As was 
noted in for a = 1, there is a critical excluded volume radius above which there is no discontinuity in the entropy 
vs. energy plot. We observed that this trend is generic for all < a < 3 and holds for both excluded volume and 
soft potential cutoffs. Now the critical cutoff radius rda) increases with increasing a, roughly varying in value from 
below 10^'^ for a = 1/4, to above I0~^ for a = 5/2, respectively. For a system with a cutoff radius larger than rJa), 
the entropy vs. energy plot is continuous and exhibits all characteristics of a normal first order phase transition [0,^: 
the convex dip and associated with it an interval of energies, where dP sie) j de^ is positive and the heat capacity is 
negative (Fig. 5). 




Fig. 5. Entropy derivative /3s(e) = ds{e)/de vs. energy e plot for a — 2 and central core radius 0.5 (solid line), and a = 1 and 
soft potential radius 0.05 (dashed line). 



We leave the more detailed study of the difference between the gravitational-like and the normal first order phase 
transitions and the nature of crossover between them for a future publications. 



VI. ACKNOWLEDGMENTS 



The authors are grateful to B. Miller and especially to H. J. de Vega for helpful discussions. This work was supported 
by the Office of Basic Engineering Science of the U.S. Department of Energy, under Grant No. DE-FG 02-88-ER 
13847. 



[1] V. A. Antonov, Vestn. Leningr Gos. Univ. 7, 13 5 (1962). 

[2] H. J. de Vega and N. Sanchez, |astro-ph/0101567| -/0101568. 

[3] T. Padmanabhan, Phys. Rep. 188, 285 (1990). 

[4] P. J. Klinko and B. N. Miller, Phys. Rev. E. 62, 5783 (2000). 

[5] V. P. Youngkins and B. N. Miller, Phys. Rev. E. 62, 5583 (2000). 

[6] M. K.-H. Kiessling, J. Stat. Phys. 55, 203 (1989). 



8 



[7] I. Ispolatov and E. G. D. Cohen, Physica A 295, 475 (2001). 
[8] R. W. Gorling and A. Hiillor, Z. Phys. B 90, 207 (1993). 
[9] D. H. E. Gross, "Microcanonical Thermodynamics", World Scientist, 2001. 
[10] F. Bavaud, Rev. Mod. Phys. 63, 129 (1991). 



9 



